library(tidyverse)
library(tidybulk)
library(ggpubr)
library(clusterProfiler)

g2r_count <- read_csv("Archive/G2R_memory_B/data/readcount.csv", name_repair = make.names)

g2r_bulkobj <- g2r_count |>
  ensembl_to_symbol(.ensembl = X) |>
  pivot_longer(
    cols = where(is.numeric),
    names_to = "sample",
    values_to = "abundance"
  ) |>
  mutate(
    group = str_extract(sample, "[a-z]+"),
    batch = str_extract(sample, "[1-3]")
  ) |>
  filter(!is.na(transcript)) |>
  select(-c(X, ref_genome)) |>
  tidybulk(sample, transcript, abundance)

g2r_scaled <- g2r_bulkobj |>
  aggregate_duplicates() |> # useful to merge transcripts in single gene, also can apply to data frame
  identify_abundant(factor_of_interest = group) |>
  scale_abundance()

# Plot count densities
g2r_scaled %>%
  ggplot(aes(abundance_scaled, group = sample, color = batch)) +
  geom_density() +
  scale_x_log10()

# Reduce data dimensionality with PCA
g2r_scaled |>
  reduce_dimensions(method = "PCA", .dims = 2, action = "get") |>
  ggplot(aes(PC1, PC2, color = group, shape = batch)) +
  geom_point(size = 5)

# Plot all-vs-all MDS dimensions
g2r_scaled |>
  reduce_dimensions(method = "MDS", .dims = 2, action = "get") |>
  ggplot(aes(Dim1, Dim2, color = group, shape = batch)) +
  geom_point(size = 5)

g2r_adj <- g2r_scaled |>
  adjust_abundance(~ group + batch) |>
  test_differential_abundance(~ group + batch,
                              method = 'edgeR_likelihood_ratio') |>
  keep_abundant() |>
  pivot_transcript() |>
  select(-c(merged_transcripts, .abundant))

g2r_adj |>
  dplyr::count(type = case_when(
    logFC > 0 & FDR < .05 ~ 'Up in RR',
    logFC < 0 & FDR < .05 ~ 'Down in RR',
    .default = 'NS'))

g2r_adj |>
  write_csv('Archive/G2R_memory_B/results/g2r_edger_LRT_res.csv')

# volcano plot -------
g2r_adj |>
  mutate(type = case_when(
    logFC > 0 & FDR < .05 ~ 'Up in RR',
    logFC < 0 & FDR < .05 ~ 'Down in RR',
    .default = 'NS'
  )) %>%
  mutate(feature = ifelse(type != 'NS', as.character(transcript), NA)) %>%
  ggplot(aes(x = logFC, y = -log10(PValue), label = feature)) +
  geom_point(aes(color = type)) +
  geom_hline(yintercept = filter(g2r_adj, FDR < .05) |> slice_max(PValue) |> mutate(pvalue = -log10(PValue)) |> pull(pvalue), linetype = 'dashed') +
  ggrepel::geom_text_repel() +
  scale_color_manual(values = c('blue','grey','red')) +
  theme_classic()

ggsave('g2r_edger_LRT_volcano.pdf',
       path = 'Archive/G2R_memory_B/figures/')

# pathway enrichment --------
g2r_up_go_ora <- g2r_sig |>
  filter(logFC > 0 & FDR < 0.05) |>
  pull(transcript) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'ALL',
           minGSSize = 3,
           readable = TRUE,
           pool = TRUE) |>
  pluck('result')

g2r_up_go_ora |>
  filter(qvalue < .05) |>
  ggplot(aes(fct_reorder(Description, Count), Count, fill = qvalue)) +
  geom_col() +
  coord_flip() +
  scale_fill_distiller(palette = 'Reds') +
  theme_pubr() +
  labs_pubr() +
  labs(x = 'Pathways')

g2r_up_down_ora <- g2r_sig |>
  filter(logFC < 0 & FDR < 0.05) |>
  pull(transcript) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'ALL',
           minGSSize = 2,
           readable = TRUE,
           pool = TRUE) |>
  pluck('result') |>
  filter((qvalue < .05 | p.adjust < .05)& Count > 1)

g2r_up_down_ora |>
  filter(qvalue < .05) |>
  ggplot(aes(fct_reorder(Description, Count), Count, fill = qvalue)) +
  geom_col() +
  coord_flip() +
  scale_fill_distiller(palette = 'Blues') +
  theme_pubr() +
  labs_pubr() +
  labs(x = 'Pathways')

g2r_gsea <- g2r_adj |>
  pull(logFC) |>
  set_names(g2r_adj$transcript) |>
  sort(decreasing = TRUE) |>
  gseGO(ont = 'ALL',
        OrgDb = 'org.Mm.eg.db',
        keyType = 'SYMBOL',
        minGSSize = 3) |>
  pluck('result')

g2r_gsea |>
  ggplot(aes(fct_reorder(Description, NES), NES, fill = qvalue)) +
  geom_col() +
  coord_flip() +
  scale_fill_distiller(palette = 'Reds') +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  labs(x = 'Pathways')

gseares <- g2r_adj |>
  pull(logFC) |>
  set_names(g2r_adj$transcript) |>
  sort(decreasing = TRUE) |>
  gseGO(ont = 'ALL',
        OrgDb = 'org.Mm.eg.db',
        keyType = 'SYMBOL',
        minGSSize = 3)

gseares |> enrichplot::gseaplot2(geneSetID = 1, title = 'GO CC:Intermediate filament')
